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Abstract 

By the use of the variational method with exponential trial functions the upper and 
lower bounds of energy are calculated for a number of non-relativistic three-body Coulomb 
and nuclear systems. The formulas for calculation of upper and lower bounds for expo- 
nential basis are given, the lower bounds for great part of systems were calculated for 
the first time. By comparison of calculations for different bases the efficiency of exponen- 
tial trial functions and their universality in respect to masses of particles and interaction 
are demonstrated. The advantage of exponential basis manifests mostly evident for the 
systems with comparable masses, though its use in one-center and two-center problems 
is justified too. For effective solution of two-center problem a carcass modification of 
the trial function is proposed. The stability of various three-particle Coulomb systems is 
analyzed. 

1 Introduction 

Among existing methods of calculation of non-relativistic bounded systems the variational 
method seems to be the most universal one as it is applied equally well for the solution of 
atomic and nuclear problems. It is essential that the variational method allows to find not 
only the upper (Eu) but also the lower (El) estimates of the energy. As to potentiality of the 
method, there are many examples of highly accurate calculations of three- and more-particle 



systems |]1|-|2T|. For instance in the three-body Coulomb problem the precision is amounted up 
to a score of decimal places. Of course in real physical systems the relativistic and other effects 
lead to corrections in the energy already in 5-7-th decimal place and therefore, in practical use 
a variational procedure which ensures a reasonable accuracy with least computational efforts 
may be acceptable. 

Historically the first variational expansion in a three-particle problem was suggested by 
Hylleraas in perimetrical coordinates in the form of exponent function multiplied by a poly- 
nomial with integer nonnegative powers. Later negative and fractional powers were added [Q], 
0, besides Frankowski and Pekeris introduced logarithmic terms. In the next this basis is 
referred as 'polynomial' one. 

Another possibility is to use a purely exponential basis. It ensures a good flexibility of 
the variational expansion due to the presence of many scale nonlinear parameters. Whereas 
the Hylleraas basis is practically oriented on solution of uni-center Coulomb problems the 
exponential basis is good for systems with any masses of particles and types of their interaction 
0. Besides, the calculations with exponential basis are more simple and uniform whereas 
for polynomial basis they became more and more complicated as number of terms increases, 
especially if the logarithmic terms are included. 
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Instead of exponents another non-polynomial functions, gaussians, can be used. They are 
not less than exponents fitted for the systems with any masses of particles, moreover they are 
applicable to systems with arbitrary number of particles. For this basis all the formulas needed 
for calculation of both the upper and lower bounds are given in paper |^, and different 3-, 4-, 
and 5-particle systems were calculated there. For the upper bound a generalization is given in 
article for arbitrary orbital moments. Nevertheless, our analysis have shown that at least 
for three-body variational calculations with not very high number of parameters the precision 
for gaussian basis is lower than for exponential one. 

Our principal goal was not the striving for improvement of existing super-high precision 
calculations but the analysis of efficiency of exponential and partly gaussian trial functions 
for evaluation of the upper and lower variational bounds. For this purpose the Coulomb and 
nuclear systems of particles with different masses and types of interaction are considered and 
the results of calculations are compared with those published in literature. 

To facilitate such a comparison we will characterize the accuracy of calculations of Eu and 
El by the values: 

Su--h(^^), *. = -lgf^i^V (1) 



Eq J \ El 

which determine the number of correct decimal places of Eu and E^ respectively, Eq being the 
exact value of the energy. 

The universality of exponential basis in respect to masses of particles allows us to analyze 
the problem of stability of different Coulomb three-particle systems. 



2 Method of calculation 

In three-particle problem it is convenient to use interparticle distances as coordinates together 
with the Euler angles describing the orientation of the triangle formed by the particles. In 
the case of central interaction the wave function of the ground state (and exited state with 
zero orbital momentum) depends only on interparticle distances, therefore the function can be 
written as: 

3 

\a) =exp(^-J2(ypRp), (2) 
p=i 

where are the nonlinear variational parameters specifying the scale of the basis function | a) , 
Ri is the distance between particles j and k where {i,j,k} is the triplet {1,2,3} or its cyclic 
permutation. In the case of gaussian basis Rp in is replaced by R^. 
It is convenient to use the notations: 

tp = cos/3p, Tp = {a\tp\b), Gp= {a\ R"^ \b), N = {a\b), 

where (3p is the angle at p-th particle in the triangle. Then simple calculations result in the 
following formula for matrix element of operator of kinetic energy, T, between states | a) and 
I b): 

{a\T\b) = Y: s;Gp - d;Tp - u'^N, 
p=i p=i 

where 

Si = ), u° = V d'} = ^/ \ 
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rrii being the mass of i-th particle and 1/M = 1/mi + l/m2 + l/rris. 

A calculation of matrix elements for the potential energy reduces to calculation of the 
integrals similar to those for kinetic energy. In particular, for the Coulomb interaction Vp = 

l/Rp. 

{a\Vp\h) = Gp. 

A calculation of lower variational estimate requires additional evaluations of matrix elements 
for operators T^, V"^ and VT. For this purpose it is convenient to introduce additional notations: 

Jpq = {a I Rp^Rg^ I b), Wpg = {a I tpR~^ \ b), Qpg = {a \ tptg \ b). 

Then, the matrix elements of operators and VpT + TVp are written as: 

3 3 

{a\T^\b) = u''u'N-J2{s;u' + sy)Gp + Y.{d;u' + dy)Tp- 

p=i p=i 

- E + sld;)Wpg + E s;slJpg + E d;d'gQpg. 



3 3 

(a I VpT + TVp I b) = E« + I I b) - Y.{dl + d\){a \ Vptg \ b) - 

q=l q=l 

- {u'' + u^){a\Vp\b). 
In the particular case of Coulomb potential: 

(a I Vp/Rg I b) = Jpg, {a \ Vptg I b) = Wgp. 

The calculation of the matrix elements of V"^ is similar to calculation of (a | | 6). In 
particular, for the Coulomb interaction a simple formula takes place: 

(a I VpVg I b) = Jpg. 

The trial function is written as a superposition of basis functions (0): 

N 

^^^ = Y.Ca\a), (3) 

a=l 

where Ca are linear parameters. 

Evidently, the difficulties arise mostly at optimization of the non-linear parameters. The 
possibilities of the deterministic procedures are soon exhausted as the number of terms in 
expansion increases. Therefore, a specially designed procedure of global stochastic searching 
was used. Briefly it is the following: at each Monte-Carlo probe a random point is chosen in 
3A^-dimensional space of non-linear parameters according to previously accepted distribution 
function. Then the coordinatewise optimization is carried on, at first the stochastic one and 
then the deterministic one. At this stage the best points are selected for subsequent detailed 
optimization. Mentioned above distribution function is found by a procedure similar to that 
described in [0. 
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3 Efficiency of calculations for various systems 



To understand better what are the possibihties of exponential basis and described above opti- 
mization procedure in calculations of systems with different masses of particles and interactions 
a number of Coulomb and nuclear systems were considered. Among them: He atom, hydrogen 
ion H~, positronium ion Ps~ (e"'"e~e~), meso-systems a^~e~, pfi~^e~ , ppfi~, /i^/i"'"e~, ^~^e~e~, 
two-center Coulomb systems ppe~ , dde~ , tte~ , as well as nuclei and 

The composition of the majority of considered Coulomb systems with the particles of unit 
charge can be presented as X^Y^Y^ , the identical particles being denoted as Y. The binding 
energies decrease together with the values of masses but the accuracy of calculation depends 
only on the ratio of masses. For these systems the upper and lower bounds were calculated with 

= 30 in expansion d), corresponding values 6u and 6l were plotted in Fig. 1 as the functions 
of mass ratio, ^yx = Tny/mx- In calculations of the lower bound the non-linear parameters 
were accepted to be equal to these found for the upper bound. 

As expected, the increase of ratio ^yx leads to the decrease of values 5u and 5l due to arising 
difficulties in description of motion of heavy particles. Nevertheless even at the approach to 
the two-center limit {^yx ^ 1) the accuracy of calculations remains still satisfactory. For 
comparison, in Fig.l the results of most detailed calculations with polynomial basis are 
presented too. It is seen that the accuracy of calculations with polynomial basis becomes 
bad for E,yx > 0.1 in spite of large values of A^. This comparison shows that the exponential 
basis is applicable for a wider range of values of ^yx than the polynomial basis. 

Note that in the case of Gaussian basis 6u and 6l decrease even more slowly than for 
exponential basis (see Fig.l) though the latter provides generally higher precision. 

The exponential basis can be used as well in the case of nuclear systems, even inconvenient 
for calculation (weakly bounded systems, short-range attractive potentials with strong repulsion 
at small distances between particles, that can be identical or not identical). As a particularly 
'inconvenient' system hypertritium, (consisting of npA), was chosen. For comparison a 
more 'convenient' three-nucleon system, ^H, was considered. In these calculations two types 
of model nuclear A^A^-potentials were used, (i) purely attractive potential A^A^ — 1 and (ii) 
attractive potential with a soft core A^A^ — 2: 



the parameters Va^r and Ra,r are given in Appendix B. 

The convergence of the upper and lower estimates for exponential basis is illustrated in 
Table 1 and in Fig. 2 for various Coulomb and nuclear systems. As seen from Fig. 2 the 
dependence of 6u and Sl on the IgA^ is close to the linear one. In accordance with Fig. 1 
the accuracy decreases as the system approaches to the adiabatic limit, and in parallel the 
convergence of variational estimates deteriorates (it is characterized by the slope of curves in 
Fig. 2). Note that the precision of calculations for considered nuclear systems is generally 
similar or even better than for Coulomb systems. 

Besides, in Fig. 2 some results of calculations with the Gaussian basis are shown (dotted 
lines). It is seen that the convergence of the upper and lower bound is similar to that for 
exponential basis whereas the accuracy is significantly lower. 



VNN{r) = Vr{Rr/R)exp{-R/Rr) - Va{Ra/R)exp{-R/Ra) 



(4) 



and attractive AA^-potential: 




(5) 
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4 Comparison of results for different bases 



Comparison of efficiency of different variational expansions is convenient to carry out on stan- 
dard systems calculated by many authors. Such systems are °°He and °°B.~ considered in 
Fig- 3 the values of 6u and 6l, are plotted for these papers where the most 



detailed calculations of atom °°He were carried out, the results of our calculations are presented 
there too. Similarly, in Fig. 4 Su and Sl are presented for hydrogen ion °°H~. 

It is necessary to emphasize that both cases are examples of one-center systems. Therefore, 
this is a reason to expect that the expansions especially designed for one-center problems 
will gain the advantage. This is generally confirmed by our analysis. Up to present the most 
accurate many-parameters calculations of °°He were carried out using polynomial or polynomial- 
logarithmic bases. As it is seen from Fig. 3 and Fig. 4 the convergence of the variational 
expansions for these bases is generally better than for exponential or gaussian bases. 

On the other hand, up to 6u ~ 12 the use of exponential basis is justified as it assures the 
same precision at lower number of terms (see Fig.3 and 4). Note that the over-high precision 
in non-relativistic calculations without taking into account relativistic and other corrections 
(that appears far before 6 ~ 12) have no physical meaning, though they are interesting from 
computational point of view. 

As to the lower bound calculations they are rare in literature and we estimate the number 
of up to which the calculations of El with exponential basis are justified (in the same sense 
as for Eu) as 100 - 200. 

Another limiting case is the adiabatic one (i.e. a two-center system with two heavy particles). 
In this case the use of polynomial basis leads to unsatisfactory results, and the exponential 
basis is evidently preferable (see Fig. 1). Moreover, the use of complex scale parameters in 
exponential basis increases significantly the accuracy of calculations [^. The most accurate 
calculations of two-center systems were carried out in the framework of the Born-Oppenheimer 
approach P3 or its modifications EM , PSI . In particular, in paper P3[ the energy of the system 



H2 {ppe ) was calculated with precision 6 ^ 12 but this is only some better than that of the 
calculation of |2^ with exponential functions (note by the way that the number of basis terms 



in was less than in ||25|| ). 

A more effective modification of exponential basis in two-center calculations is: 

\a) = exp(-a^i?i - a^i?2 - "3^3 - P^RI), (6) 

where is a distance between the heavy particles. Note that the dependence of this function 
on R3 can be presented as exp(— /5°(i?3 — R^Y), where R^ is the new variational parameter 
connected with ag. Note that basis (P) is, in a certain sense, a particular case of 'carcass' 
functions (constructed on the base of gaussians in paper ^^), whose use together with gaussians 



might be useful in nuclear physics for calculation with potentials changing the sign. 

For functions (^) all the integrals needed for calculations of the upper variational bound are 
expressed in the closed form in terms of conventional functions. For instance, the basic integral 
can be calculated as: 



1 ^-^,R^-mldR^ Je-^^^dR2 f e-^^^^dR, = , , ^ , (7) 

J J J y P ^1 — ^2 



R2+R3 

qOOO — '' 

\R2-R3\ 

rOO — /2 



where F{z) = e dt. 

The calculations of the ground state of the system ppe~ with this modified basis lead to 
significantly better results than with purely exponential or gaussian bases. In particular, in our 
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calculations it has been shown that even a single function (j^) provides a better precision than 
50 exponents or gaussians. Moreover, the basis (0) is more flexible than the exponential basis 



with complex parameters used in |2^. For instance, the result of calculations with = 20 
for ppe~ turns out to be better than that of paper with 200 complex exponents (1400 
variational parameters) and better than calculations with 300 functions for systems fi^fi^e~, 
dde~ and tte~ . 

In addition to the preceding discussion of two limiting cases (one- and two-center problem) 
it is necessary to indicate that there exists a large region of values of ^yx between 10~^ and 
10^ where the exponential basis is beyond compare. Note that this is the region where the 
great part of known three-particle Coulomb systems is located. Thus, apart from gaussians, 
the exponential basis seems to be the most universal one in comparison with other approaches, 
applicable equally well to Coulomb and nuclear three-particle systems. 



5 Stability of Coulomb Systems 

All considered above Coulomb systems except two (ae^e^ and afj,~e~) had summary charge 
±1 and consisted of three single- charged particles from which two are identical. All systems of 
such type are stable in respect to separation of one of the particles. However this is not the 
case for other type of three particle Coulomb systems. For analysis of stability of Coulomb 
systems and for calculation of their energy it is natural to use the variational procedure with 
exponential basis as it is most universal in respect to masses of particles (see also [p7[| ). 

In general case the structure of a Coulomb system of three single-charged particles with total 
charge ±1 may be presented in the form X^Y^Z^ where my < m^. The stability of the system 
depends on two ratios of masses, ^yx = mY/mx and C,zx = ^z/^x- A boundary delimiting 
the regions of stable and unstable systems is determined from the condition of coincidence 
of the energy of the three-particle system with that of the two-particle system X'^Z'^ . The 
corresponding equation determining the interdependence between ^yx and ^zx can be written 
as: 

JKKyx.Kzx) = ^^x^Z^) 

The solution of this equation is presented in Fig. 5 by the curve A. It is seen that not only 
systems with two identical particles are stable but also two-center systems (two heavy particles 
with identical charges plus light particle with opposite charge). In contrast, a system containing 
two heavy particles of opposite charges are unstable. An exception can occur if all three particles 
have nearly equal masses. This takes place for instance for exotic systems p'^p'S^ {C,yx = 1? 
^zx = 1.2749), /i+/i-7r± (^yz = 1, ^zx = 1.3213) and 7r+7r-/i± (^yx = 0.7568, ^zx = 1) for 
which / = 0.008745, 0.006069 and 0.002354, respectively. Of course, a three-particle system 
which is stable with respect to emission of one of the constituent particles can be unstable 
in the excited state. This problem was considered in ||^ for symmetric (XYY) systems with 
mY/mx <^ 1. 

For the case of systems of the type x~^"^Y^'^Z~'^ containing multiple-charged particles the 
situation is quit similar to the case of single-charged particles considered above. Among three- 
body systems containing single and double charged particles the systems of the type X^^Y~ Z~^ 
and X~^^Y^ Z^^ are unstable at any ratio of their masses, whereas the systems X^^Y^ Z^ are 
always stable. As to the systems of the type X^^Y Z^ they can be stable only for restricted 
values of ratios of their masses. The corresponding boundary is shown in the same Fig. 5, 
curve B. 
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A Standard integrals 

A calculation of matrix elements of the Hamiltonian and its square reduces to the evaluation 
of the following integrals: 



The integrals I^^"^{xi,X2,X3) with non-negative indexes are the uniform polynomials of the 
{k + I + m + 3)-th degree with respect to the variables Ai = 1 / {xi + X2 + — Xi). 
To calculate the upper variational estimate the following integrals are necessary: 



N = 2/000((Ai + A2)(A2 + A3)(A3 + Ai)-AiA2A3); 

Gi = I^^'^{A,A2 + A2A, + AsA^ + 2Aly, 

Ti = /1"-4J00M2A3(A2 + A3). 



(Here and further an unimportant numerical factor IGvr^ is dropped.) 

For presentation of the integrals (|^) with negative indexes it is convenient to use the following 
notations: 



Bi = {x2-X3) ^ B2 = {x3-Xi) \ B3 = {xi-X2) 



To calculate the lower variational estimate the following integrals are necessary: 



Jl2 = J21 = /°°°(Ai+A2); 

Jii = A?5™-Ai5p^+/ooo(Ai/2 + (Ai + A2 + A3)/2); 

W^2 = /°°Mi(A3 + 2^2)+ 2^35^-2^24=^!; 

1^21 = /°°M2(A3 + 2A^) - 2AIS?^ - 2AlSf^; 

M/n = G1-2/00M2A3; 

Qi2 = Q2i = iV-4/ooM3K + A2)-8A35P; 

Qii = 3^2 — 25*^2^2 + 35*^2^2 ~ 35*^2^2 + 25^2) + 

+ 3^3 ( - 24V3 - 35a^i + 34V3 + 241) + 
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r5-i-i 



60<^ xr^[iVf - ^2/13(^3/4 + ^2^3/6 + ^2^2/6 + All A 



+ 



+ Xi 



+ Xi 



A^f - ^2^3(^2/2 + A3/2)] + xi 



+ 



+ xr^[Li2(l-f) + Li2(l-^ 



2 A3 ^ 6 



The expression for 1 1 contains the di-logarithmic function Li2(2) = — ^^t ^ ln(l — t)dt. 
If u = X1/X2 and V = Xi/xs are simultaneously small one can use for it the expansion: 



r5-l-l 



w 



(n + 5)! 



(^ + 5)! 



P,(«) = Pg_2(a) 



?i=0 



n!(n + 6) n\[n + 6) 



q = 2,3, 



Po = -ln(l + a), Pi = l, w = -. 

u 



These formulas are used if max(u,v) < 0.3. 



B Model Nuclear potentials 

Parameters of used nuclear model potentials and (|^) are given in Table 4. 

In calculations of hypertritium A^A^-potentials NN-1 and NN-2 were used. The radial 
parameter of purely attractive potential NN-1 was chosen corresponding to one-pion exchange 



whereas for NN-2 potential the values of Rr and Ra were adopted from paper The 
depth parameters for potentials NN-1 and NN-2 were matched to correct deuteron energy, 
additional experimental data in fitting of parameters for NN-2 potential were deuteron radius 
and phases of S-wave triplet np-scattering up to energy 300 MeV. In calculations of tritium 
the potential NN-3 was used with the same radial parameter as for potential NN-1, whereas 
depth parameters was chosen to describe the correct tritium binding energy in calculations with 
= 100 in expansion (j^). 

The radius of AA^-potential was adopted from paper while the depth parameter provided 
the correct hypertritium binding energy (Ba = 0.13 MeV) in calculations with = 100. 
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Figure 1: Dependence of 6u and on mass ratio for Coulomb systems X^Y Y 
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System 


N 


Eu, au 


El, au 


Comment 




10 


-2.903 723 6 


-2.903 83 






30 


-2.903 724 373 


-2.903 725 8 






50 


-2.903 724 375 9 


-2.903 725 2 






100 


-2.903 724 377 009 


-2.903 724 414 






200 


-2.903 724 377 030 3 


-2.903 724 391 






300 


-2.903 724 377 033 2 


-2.903 724 380 




aee 


300 


-2.903 304 557 732 3 


-2.903 304 561 


m« = 7294.2996me 


°°H- 


10 


-0.527 750 546 


-0.528 062 






30 


-0.527 751 009 425 


-0.527 764 






50 


-0.527 751 015 895 


-0.527 752 977 






100 


-0.527 751 016 400 


-0.527 751 663 




pee 


100 


-0.527 445 880 971 


-0.527 446 533 


nip = 1836.1527me 


ixee 


100 


-0.525 054 806 098 


-0.525 055 501 


m^, = 206.768262me 


Ps- 


10 


-0.262 003 563 


-0.262 744 






30 


-0.262 005 053 


-0.262 026 






50 


-0.262 005 068 6 


-0.262 008 7 






10 


-0.494 374 


-0.495 7 


In mcso-atomic units 




30 


-0.494 386 645 


-0.494 408 


rup = 8.8802444m/, 




50 


-0.494 386 790 


-0.494 391 1 


ddix 


10 


-0.531 044 


-0.534 4 


In meso-atomic units 




30 


-0.531 109 463 


-0.531 241 


nid — 17.7516751m^ 


tt/i 


10 


-0.546 224 


-0.551 2 


In meso-atomic units 




30 


-0.546 371 871 


-0.546 517 


mt = 26.5849388m^ 


lUie 


10 


-0.583 276 


-0.604 5 






30 


-0.584 757 


-0.588 82 






50 


-0.584 995 


-0.586 267 






20 


-0.585 126 081 8 




'Carcass' basis 


ppe 


10 


-0.591 03 


-0.625 






30 


-0.595 02 


-0.612 






50 


-0.595 67 


-0.606 






20 


-0.597 139 058 5 




'Carcass' basis 


dde 


10 


-0.591 38 


-0.621 






30 


-0.596 06 


-0.610 






20 


-0.598 788 780 3 




'Carcass' basis 


tte 


10 


-0.591 59 


-0.615 






30 


-0.596 34 


-0.608 






20 


-0.599 506 906 3 




'Carcass' basis 


PHe 


10 


-0.584 18 


-0.645 




a/ie 


10 


-1.947 287 542 


-1.947 429 


In meso-atomic units 




30 


-1.947 287 553 22 


-1.947 290 320 


Eo =-1.947 287 553 40, 
rria = 35.2776559m^, 
me = 0.00483633218m^ 


Q T T 

aH 


50 


-2.359 478 5 


-2.437 6 


Potentials NN-1 and AN-1; energies 

are in MeV 


iH 


50 


-2.358 597 8 


-2.808 


Potentials NN-2 and AN-2; energies 
are in MeV 




50 


-8.480 037 312 


-8.480 045 6 


Potential NN-3; energies are in MeV 



Table 1: Upper and lower bound|^or Coulomb and nuclear systems 



Variant 


K, MeV 


Rr, Fm 


Va, MeV 


Ra, Fm 


NN-1 








50.6414 


1.4 


NN-2 


2719.20 


0.32 


730.24 


0.65 


NN-3 








40.0419 


1.4 


AN-1 








687.00 


0.23 


AN-2 








711.00 


0.23 



Table 2: Parameters of nuclear model potentials 
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Figure 2: Dependence of Su and Sl on number of terms in variational expansion for exponential 
basis 
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Figure 3: Su and 6l in calculations of atom °°He with different variational expansions. Markers 
refer to the first author of the corresponding paper. 
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